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Oxygen  reduction  at  the  mixed  ionic  and  electronic  conductive  (MIEC)  SOFC  cathodes  occur  along  a 
surface  pathway  with  oxygen  transport  through  a  triple-phase  boundary  (3PB),  and/or  a  bulk  path¬ 
way  via  bulk  cathode  and  electrolyte/cathode  interface  (2PB).  The  identification  of  the  rate-limiting 
steps  (RLSs)  for  each  path  and  understand  on  their  interactions  are  important  to  the  SOFC  cathode 
kinetics.  In  this  research  a  modified  one-dimensional  continuum  model  is  developed  to  analyze  the 
oxygen  reduction  on  LSM-type  MIEC  cathode  by  incorporating  multi-step  charge-transfer  into  the  bi¬ 
pathway  kinetics.  Finite  control-volume  method  is  used  to  simulate  the  pathway  kinetic  competition, 
and  a  parametric  study  is  performed  with  different  equilibrium  concentrations  of  surface  oxygen  ion 
(Qr.eq)  and  hulk  oxygen  vacancy  (Cv,MiEc,eq).  The  I-V  profiles  show  kinetic  transition  from  3PB-control 
to  2PB-control  at  overpotentials  from  -0.2  V  to  -0.4  V,  and  the  active  reaction  zone  in  concentration 
profiles  expands  from  3PB/2PB  interfaces  by  2-4  (Jim.  The  exchange  and  local  3PB  currents  recognized 
from  simulation  indicate  limit  of  surface  oxygen  diffusion  as  the  mechanistic  process  for  pathway 
transition.  The  results  are  compared  to  the  reported  experimental  findings  to  examine  the  model’s 
assumption  of  oxygen  reduction  scenario,  and  its  potential  implication  on  SOFC  cathode  performance 
improvement. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

For  decades,  solid  oxide  fuel  cells  (SOFCs)  have  been  widely 
investigated  as  a  promising  energy  conversion  system  with 
advantageous  efficiency,  fuel  flexibility  and  emissions  [1,2].  Yittria- 
stablized-zirconia  (YSZ)  is  commonly  applied  as  a  stable  and 
economical  electrolyte  material  [3,4],  while  various  perovskite 
oxides  with  mixed  ionic  and  electronic  conductivity  (MIEC)  are 
developed  into  cathode  materials.  Among  them  strontium-doped 
lanthanum  manganite  (LSM)  appears  as  the  conventional  choice 
due  to  good  electrical  property,  reliability  and  low-cost  [5]. 
It  also  serves  as  a  stand-point  for  understanding  of  the  reac¬ 
tion  fundamentals  and  development  of  more  active  cathodes, 
such  as  (La,Sr)Co03  (LSC)  and  (La,Sr)(Co,Fe)03  (LSCF).  It  is  eco¬ 
nomically  desirable  for  the  SOFC  to  operate  at  an  intermediate 
temperature  typical  of  700-850  °C,  but  higher  electrode  per¬ 
formance  is  then  required  to  overcome  the  increased  ohmic 
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resistance.  As  a  major  portion  of  polarization  loss  is  from 
relatively  sluggish  oxygen  reduction  [6,7],  intensive  efforts 
have  been  devoted  to  investigate  and  improve  SOFC  cathodes 
kinetics.  Even  so,  questions  remain  on  the  assignment  of  rate- 
limiting  steps  (RLSs)  to  the  available  kinetic  pathways  and 
their  expected  connection  to  cathode  activation  [8,9].  Particularly 
for  the  LSM-based  cathode,  oxygen  transport  across  the  cath¬ 
ode/electrolyte  interface  is  known  to  be  restricted  by  the  low 
intrinsic  ionic  conductivity  of  LSM.  The  reported  RLSs  accord¬ 
ingly  include  charge-transfer  to  the  oxygen  adsorbate,  surface 
oxygen  diffusion,  and  incorporation  of  oxygen  adsorbate  at  the 
triple-phase-boundary  (3PB)  as  suggested  by  AC  impedance  and 
DC  polarization  study  [10-13].  When  conditioned  by  cathodic 
current,  partial  reduction  of  Mn  cations  and  subsequent  oxygen 
vacancy  formation  in  LSM  are  often  suggested  as  the  activa¬ 
tion  mechanism  [14-16],  although  removal  of  surface  passive 
species  and  microstructural  change  are  also  considered  respon¬ 
sible  [17-19].  Nevertheless,  debates  between  chemical  versus 
electrochemical  control  on  the  gas-electrode  reactions  further 
complicate  identification  of  RLSs  and  interpretation  of  activation 
behaviors  [20,21  ].  In  a  sense  of  reducing  experimental  uncertainty 


M.  Gong  et  al.  /  Journal  of  Power  Sources  201  (2012)  204-218 


205 


Nomenclature 

C0-  &  C0-  concentration  of  adsorbed  surface  oxygen  ion 

ad 

C0-  ea  &  C0-  equilibrium  concentration  of  adsorbed  sur- 

ad’  4  ’  4 

face  oxygen  ion 

C0-  3pB  concentration  of  adsorbed  surface  oxygen  ion  at  3PB 

zd’ J 

(x  =  0) 

Cv.miec  &  Cv  concentration  of  oxygen  vacancy  inside  MIEC 
cathode 

Cv.MiEc.eq  &  Cv,eq  eQuilibrium  concentration  of  oxygen 
vacancy  in  cathode 

CV)miec,x=1c  concentration  of  oxygen  vacancy  at  MIEC/EC 
boundary  (x  =  lc) 

Cv.ysz  concentration  of  oxygen  vacancy  inside  YSZ  elec¬ 
trolyte 

£eqE  working  electrode  potential  at  thermal  equilibrium 
(OCV) 

Ere  reference  electrode  potential 

Eocv  open-circuit  potential  for  the  cell 

E2 pb  2PB  electrode  potential  as  electrostatic  potential 

difference  across  2PB 
2PB  electrode  potential  at  OCV 
E3PB  3PB  electrode  potential  as  electrostatic  potential 
difference  across  3PB 

Es  surface  electrode  potential  as  electrostatic  potential 

difference  across  gas/MIEC 
0MIEC  bulk  electrostatic  potential  of  MIEC 

0S.miec  surface  electrostatic  potential  of  MIEC 

0YSZ  bulk  electrostatic  potential  of  YSZ 

^miec  electrochemical  potential  of  species  i  in  MIEC 
^fchem  chemical  potential  of  species  i  in  MIEC 
fiYSZ  electrochemical  potential  of  species  i  in  YSZ 
/zfdiem  chemical  potential  of  species  i  in  YSZ 


and  fundamental  ambiguity,  mechanistic  approaches  merit,  either 
by  using  patterned  micro-electrodes  to  observe  RLSs  as  func¬ 
tions  of  well-defined  electrode  geometric  parameters  [22-24]  or 
through  numerical  parametric  simulations  to  study  the  fundamen¬ 
tal  kinetics  of  the  reaction  pathways  and  facilitate  experimental 
design. 

Presently  there  are  two  different  pathways  recognized  for 
oxygen  reduction  on  an  MIEC  cathode  [25],  in  which  overall 
oxygen  transport  occurs  through  exposed  3PB  or  the  buried 
cathode/electrolyte  boundary  (2PB).  Concurrent  mass  and  charge 
transfer  proceed  as:  (i)  gaseous  diffusion  to  the  electrode  surface; 
(ii)  surface  oxygen  exchange  reactions;  (iii)  surface  and  bulk  solid- 
state  diffusion  of  oxygen  species;  and  (iv)  charge-transfer  at  the 
3PB  and  2PB.  Fundamental  models  of  the  MIEC  cathodes  have 
systematically  studied  kinetic  contribution  from  reaction  steps 
(i)-(iv)  to  unveil  the  essential  RLSs  for  oxygen  reduction  via  dif¬ 
ferent  mechanisms,  however  only  one  single  kinetic  pathway  is 
usually  considered  dominant  in  practice.  For  study  focusing  on 
2PB-pathway,  a  continuum  model  was  first  proposed  by  Alder, 
Lane  and  Steele  to  analyze  oxygen  reduction  on  single-phase  MIEC 
cathode  [20],  and  the  model  based  on  chemical  kinetics  accu¬ 
rately  interpreted  impedances  for  highly  ionic  conductors,  but  was 
unconvincing  for  poor  ionic  conductors  on  which  3PB-pathway 
reaction  may  dominate.  Sogaard  then  modeled  the  2PB-pathway 
reactions  for  a  composite  SOFC  cathode  from  similar  chemical 
kinetics  assumption  with  a  finite  volume  method  [26].  In  contrast, 
Fleig  derived  a  model  for  2PB-pathway  kinetics  by  assum¬ 
ing  electrochemical  driving  force  of  the  surface  potential  step 


from  a  dipole  layer  of  surface  oxygen  ion  [27],  and  Mebane 
further  developed  Fleig’s  approach  into  a  2-D  model  from 
first  principle  analysis,  which  not  only  compared  the  impor¬ 
tance  of  the  electrochemical  driving  force  versus  the  chemical 
one  for  2PB-pathway  reactions,  but  also  identified  geometric 
charge  inhomogeneity  [28].  On  the  other  hand  for  3PB-pathway, 
Tanner  and  Virkar  built  up  both  1-D  and  2-D  models  to  quan¬ 
tify  microstructural  effects  on  the  total  3PB  reaction  kinetics 
[29,30]. 

The  results  together  promote  understanding  and  design  of 
the  SOFC  cathode;  however  there  are  frequent  cases  that 
two  kinetic  pathways  may  co-exist  competitively,  especially 
when  a  poor-MIEC  cathode  is  can  be  transformed  into  MIEC 
by  operation  at  high  overpotential.  Takeda  and  Siebert  first 
used  semi-empirical  models  to  summarize  kinetic  contributions 
from  both  pathways  for  activation  of  LSM  cathode  [14,16]. 
Liu  then  developed  a  comprehensive  reaction  model  from  the 
bi-pathway  kinetics  [31,32].  Svensson  proposed  a  1-D  model 
to  parametrically  simulate  the  bi-pathway  kinetics  on  MIEC 
cathodes  [33,34].  Coffey  later  improved  the  model  to  estab¬ 
lish  the  3PB-to-2PB  pathway  transition  by  quantifying  the 
current-overpotential  response  [35].  These  efforts  employ  mea¬ 
surable  or  tabulated  material  property  parameters,  and  thus 
guide  experimental  probe.  However,  if  applying  the  models  on 
a  poor-MIEC  cathode  for  which  pathway  transition  most  likely 
to  appear  during  activation,  two  challenges  may  occur.  First, 
the  assumptions  on  chemical  driving  force  for  the  oxygen-MIEC 
surface  reactions  and  single-step  charge-transfer  for  interfa¬ 
cial  reactions  need  reconsideration.  Experimental  study  by  van 
Heuveln  and  other  researchers  indicated  charge-transfer  on 
LSM-type  cathode  was  a  multiple-step  process  [10-12].  Phe¬ 
nomenological  and  first-principle  analysis  also  suggested  surface 
reactions  to  be  electrochemically  driven  [27,28,36].  Second  as 
Adler  stated  [5],  simulation  of  steady-state  electrode  behavior  can¬ 
not  be  used  to  predict  impedance  and  other  transient  physical 
processes. 

Inspired  by  prior  works,  a  one-dimensional  diffusion  model 
is  developed  in  this  study  to  simulate  the  bi-pathway  oxygen 
reduction  on  a  MIEC  cathode  from  multi-step  charge-transfer 
kinetics.  This  mode  would  complement  Svensson  and  Coffey’s 
models  in  the  above  two  aspects  mentioned:  first  we  treat 
the  overall  oxygen  reduction  as  a  multiple-step  charge-transfer 
process,  in  which  not  only  the  2PB  and  3PB  interfacial  reac¬ 
tion  rates,  but  also  the  gas/MIEC  surface  reaction  rates  are 
phenomenologically  applied  with  Butler-Volmer  type  expres¬ 
sions  using  electrochemical  driving  forces.  Second,  this  model 
analyzes  the  unsteady  diffusion  equations  for  MIEC  cathode 
and  utilizes  a  time-discretization  solver  to  reach  steady-state. 
This  approach  examines  the  kinetic  pathway  transition  from 
a  more  realistic  scenario  involving  multi-step  charge  trans¬ 
fer  and  kinetic  evolution.  A  finite  control-volume  method  is 
used  to  visualize  the  transient  defect  profiles  in  the  electrode 
as  functions  of  overpotentials  and  material  physical  parame¬ 
ters.  The  assumption  of  uniform  Fermi-level  and  perfect  current 
collection  throughout  the  cathode  might  be  problematic  for 
the  systems  with  high  oxygen  vacancy  concentration  and  cer¬ 
tain  thin-film  structure  [27,37],  therefore  this  simulation  is 
primarily  designed  for  a  LSM-type  poor-MIEC  during  activa¬ 
tion.  As  the  model  focuses  on  comparison  of  oxygen  reduction 
pathways  under  charge-transfer  RLSs,  microstructural  diffusion 
effects  are  less  addressed.  The  simulation  results  reveal  that 
the  nature  of  the  3PB-to-2PB  kinetics  transition  involves  oxy¬ 
gen  vacancy  evolution  and  diffusion  limit  for  surface  oxygen 
as  2PB  kinetic  contribution  increases.  Effects  of  surface  and 
bulk  material  parameters  on  kinetic  pathway  transition  are  also 
discussed. 
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Fig.  1.  Three  mechanisms  for  oxygen  reduction  on  SOFC  cathode. 


2.  Physical  model 

2. t.  General  schemes  for  oxygen  reduction 

Fig.  1  displays  the  oxygen  reduction  on  a  SOFC  cathode  divided 
into  three  mechanisms.  Mechanism  I  depicts  materials  with  poor 
ionic  conductivity  and  essentially  a  metallic  electronic  conduc¬ 
tor  (EC),  which  represent  LSM-type  cathode  operated  under  low 
overpotential  prior  to  activation.  The  only  active  pathway  is  3PB 
charge-transfer.  Mechanism  III  describes  a  MIEC  with  good  ionic 
conductivity,  which  may  correspond  to  LSCF-type  cathode  with 
incorporation  and  diffusion  of  oxygen  adsorbate  primarily  through 
the  2PB  bulk  pathway.  The  present  study  interest  is  mechanism 
II  where  electrode  kinetics  proceeds  on  two  parallel  paths.  Such 
a  mechanism  applies  to  LSM-type  cathode  after  activation,  when 
the  portion  of  cathode  bulk  near  3PB/2PB  is  transformed  from  EC 
to  MIEC  by  oxygen  vacancy  exchange.  Oxygen  reduction  via  2PB 
pathway  then  becomes  active  and  complements  the  existing  3PB 
pathway  kinetics. 

2.2.  Multi-step  charge  transfer  model 

Oxygen  transport  on  an  LSM-type  cathode  in  contact  with  a 
YSZ  electrolyte  is  described  by  a  1-D  physical  model  in  Fig.  2. 
The  electrode  consists  of  two  layers:  (1)  the  MIEC  layer  (acti¬ 
vated  portion  of  LSM  with  ionic  conductivity)  next  to  the  3PB/2PB 
interface,  on  which  reaction  steps  (S1-S4,  B3-B4)  of  oxygen 
adsorption,  reduction,  and  diffusion  occur  via  both  3PB  and 
2PB  paths;  (2)  an  EC  outlayer  (inactivated  portion  of  LSM)  for 


surface  oxygen  adsorption  (SI).  The  boundary  separating  the 
two  layers  is  defined  by  where  the  bulk  flux  C/difr.v)  of  oxygen 
vacancy  in  the  MIEC  (Vq  MIEC)  generated  via  2PB  exchange  reaction 
approaches  zero.  The  only  geometric  factor  is  the  volume  spe¬ 
cific  surface  area  of  electrode,  which  is  denoted  as  AS/AV  in  the 
model. 

The  surface  (3PB)  pathway  consists  of  four  elementary  steps 
denoted  sequentially  as  S1-S4,  with  corresponding  physical  pro¬ 
cesses  listed  in  Table  1.  The  bulk  (2PB)  pathway  shares  the 
first  two  reaction  steps  (SI  and  S2)  with  the  3PB  pathway,  fol¬ 
lowed  by  two  additional  steps  of  B3  and  B4.  Adopting  Fleig  and 
Mebane’s  definition  [27,28],  reaction  B3  is  derived  by  break¬ 
ing  the  overall  surface  oxygen/MIEC  exchange  reaction  Oad  + 
Vq  miec  +  2e  Oo(MIEC  into  two  half  electrochemical  reactions 
S2  and  B3.  B4  is  the  net  charge-transfer  reaction  for  the  2PB 
pathway.  Above  reaction  steps  are  generally  set  from  individual 
breakdown  of  the  overall  3PB/2PB  oxygen  reductions,  by  follow¬ 
ing  previous  study  on  MIEC  cathode  kinetics  [10-13,16,27,28]. 
Although  other  kinetic-step  scenarios  may  also  merit,  the  scope  of 
this  paper  rather  tries  to  provide  a  numerical  treatment  method 
to  analyze  multi-step  charge  transfer  along  different  kinetics 
pathways. 

The  derivation  of  reaction  rates  for  each  elementary  reaction 
step  follows  the  detailed  procedures  adopted  by  van  Eleuveln, 
Mebane  and  Bockris  according  to  the  multi-step  electrode  kinet¬ 
ics  and  transition-state  theory  (TST)  [10,28,38,39].  The  kinetic 
analysis  is  based  on  the  fundamentals  that  for  charge-transfer 
reactions  multiple  electrons  (holes)  cannot  be  transferred  simul¬ 
taneously  due  to  prohibitive  activation  energy  barrier  [40]. 


lc  (total  length  simulated) 


LCatliode  (total  cathode  thickness) 


Fig.  2.  Schematic  description  of  bi-pathway  dominated  oxygen  reduction  on  MIEC  cathode  (LSM). 
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Table  1 

Electrode  processes  assigned  to  each  elementray  steps  in  surface  and  bulk  pathways. 
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Process  description 

Dissociate  oxygen  adsorption  (S:  surface  adsorption  site) 

Oxygen  intermediate  formation  (1st  step  charge-transfer) 

Surface  diffusion  towards  3PB 

3PB  incorporation  of  oxygen  intermediates  (2nd  step  charge-transfer) 
Reaction  of  vacancy  from  MIEC  bulk  with  surface  oxygen  intermediates 
2PB  exchange  between  YSZ  and  MIEC(S:  surface  adsorption  site) 


Generally  for  a  thermal-activated,  single-step  charge-transfer  reac¬ 
tion  OX  +  e*>  RE,  the  forward  rate  constant  /<R  and  backward  rate 
constant  k0  are  modified  by  electrostatic  potential  difference  Eox/re 
across  reaction  interface  as: 

kR  =  k°  exp  (-^Eox/re)  (1) 

ko  =  exp  (— E0x/re)  (2) 

where  Fis  the  faradic  constant,  a  the  symmetry  coefficient,  R  the  gas 
constant  and  T  the  reaction  temperature;  and  are  collective 
constant  sets  with  transmission  coefficient,  Boltzman’s  constant, 
standard  chemical  potentials  and  activity  coefficients  as  derived 
by  Mebane  and  Liu  [28]. 

Following  this  formulism  the  net  rates  of  the  consecutive  reac¬ 
tion  steps  for  surface  (3PB)  and  bulk  (2PB)  pathways  can  be 
phenomenologically  expressed  with  species  concentrations  as: 


rsi  =  ksiPo^r(1  -9)-  ID  (3) 

rS2  =  kS2re  expi-oisfEs)  -  k~2C0-  exp[(l  -  a s)fEs]  (4) 

rS3  =  l<S3 C0-  -  k~3C0-  3pB  (5) 

ad  zd 

rS4  =  kS4C0-  3PBCV)ysz  exp(-Q'3pB/E3pB)  -  k~4r{  1  -  6) 

zd’ 

x  exp[(l  -  «3pb)/E3pb]  (6) 


fB3  =  fcB3 c0-  Cv.miec  exp(as/Es)  -  fcB3r(l  -  9)  exp[-(l  -as)fEs\ 

ad 

(7) 

rB4  =  /<b4Cv,ysz  exp(-2a2PR/E2PB)  -  kB4 Cv.miec 

x  exp[2(l  -  o'2pb1/E2pb ]  (8) 

where  ks,  k~,  kB  and  k~  are  the  forward  and  backward  rate 
constants  for  the  denoted  reaction  steps  in  the  two  pathways. 
C0-  ,  Cy, miec  and  CV)ysz  are  the  local  concentrations  for  sur- 

ad 

face  oxygen  ion,  oxygen  vacancy  inside  MIEC  cathode  and  oxygen 
vacancy  of  the  electrolyte,  respectively;  r  is  the  surface  adsorp¬ 
tion  site  density,  and  0  is  the  fraction  of  oxygen  site  coverage.  The 
symmetry  factors,  a3PB,  a2PB  and  as  are  taken  as  0.5,  and  thermal 
factor /is  given  by /=F/RT.  Note  that  although  Eq.  (5)  is  written  in 
a  way  resembling  a  kinetic  reaction  for  formulism  convenience, 
the  process  of  step  S3  is  actually  surface  oxygen  diffusion.  And 
the  rate  constants  in  Eq.  (5)  are  the  same  as  a  ratio  of  the  sur¬ 
face  diffusivity  over  effective  surface  diffusion  distance  according 


to  Fields  1st  law  [10].  Elementary  steps  S2,  S4  and  B3  are  simple 
charge-transfer  reactions,  while  step  B4  in  Eq.  (8)  is  also  a  single- 
step  charge-transfer  reaction.  Following  Mebane’s  practice  [28], 
the  electric  potential  terms  E  above  are  the  electrostatic  poten¬ 
tial  differences  across  respective  interfaces,  with  Es  as  the  surface 
electrode  potential,  F3PB  the  3PB  electrode  potential  and  F2pB  the 
2PB  electrode  potential.  Detailed  definitions  and  assumptions  of 
these  electrode  potentials  are  separately  discussed  below,  and  can 
be  further  found  in  Appendix  B. 

2.3.  Potential  and  overpotential  (net  driving  force)  for  the 
reactions 

The  potential  and  overpotential  of  the  MIEC  electrode  is  con¬ 
sidered  in  an  electrochemical  cell  with  a  3-electrode  configuration 
with  the  cathode  as  working  electrode  (WE),  a  Pt  counter  electrode 
(CE)  and  a  reference  electrode  (RE).  At  open  circuit,  the  net  reaction 
rate  is  zero  for  the  WE.  This  yields  the  Nernst  relationship  between 
the  equilibrium  potential  of  WE  and  oxygen  partial  pressure  no 
matter  which  step  is  presumably  the  RLS  [10]: 

E^e  =  constant  +  ^  In  p^E  (9) 

The  open  circuit  potential  of  the  cell  (F0Cv)  is  the  potential  differ¬ 
ence  between  the  WE  and  RE,  and  will  only  depend  on  the  oxygen 
activity  difference  according  to  Nernst  equation  [35]. 

Eocv  =  E™  -  Ere  =  -  g  in  pEE  +  g  in  p£E  (10) 

Evaluated  against  an  ideal  RE,  the  measured  overpotential  is  due 
to  the  change  in  WE  potential: 

Ptot  =  E-£0cv  =  EWE-EWE  (11) 

Next  individual  electrode  potential  and  overpotential  at  3PB, 
2PB  and  cathode  surface  need  to  be  evaluated.  Here  distribution 
of  lateral  potential  and  Fermi  level  is  assumed  to  be  uniform  for 
MIEC  cathode,  which  is  also  taken  as  reasonable  by  other  1-D 
models  [34,35].  With  this  assumption  and  from  the  standpoint  of 
good  electronic  conductivity,  surface  electrode  potential  Es  and  3PB 
electrode  potential  F3PB  are  then  considered  equal  to  F2PB  in  Eqs. 
(4)  and  (6)-(8).  In  this  model  although  we  follow  the  Svensson 
and  Coffey’s  assumptions  on  electrode  potentials  in  that  F3PB  =  F2PB 
[33,35],  we  rather  adopt  Mebane’s  method  of  TST  kinetic  analysis 
to  define  F3PB  and  E2PB  [28],  which  are  the  interfacial  electrostatic 
potential  differences,  as  the  driving  forces  for  the  charge-transfer 
reactions,  because  this  difference  represents  activation  energy 
actually  driving  the  reaction.  Hence  the  net  driving  forces  for  3PB 
and  2PB  reactions  can  be  defined  as  the  measured  the  electrode 
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overpotentials  from  OCV,  while  omitting  anode  overpotential  and 
cell  ohmic  potential  loss: 

*73PB  =  ^3PB  -  ^OCV  (12) 

11 2PB  =  ^2PB  -  ^OCV  (13) 

which  are  different  from  Svensson  and  Coffey’s  models  as  the 
change  of  interfacial  electrostatic  potential  differences.  Specific 
derivation  of  r]3pp  and  ^2pb  as  net  driving  forces  from  model 
assumptions  is  shown  in  Appendix  B.  With  E3pb  =  E2pb  assumed 
for  a  1-D  model,  ^3pb  =  ^72pb  can  be  then  derived  and  the  two 
pathway  reactions  can  have  the  same  net  driving  force.  Simi¬ 
larly,  we  observe  that  ijs  =  r]3pp  =  t]2pb •  It  needs  to  be  pointed  out 
that  the  surface  potential  change  rjs  due  to  charge  build-up  from 
oxygen  surface  exchange  can  be  rather  different  from  the  cath¬ 
ode/electrolyte  interfacial  overpotential,  and  this  could  bring  about 
a  complicated  dependence  of  qs  on  r\2 pb  [27].  This  will  be  fur¬ 
ther  discussed  in  Section  5,  and  the  present  model  assumes  a 
simplified  relationship  between  rjs  on  r]2pB  for  computation  conve¬ 
nience.  Further  discussion  on  relationship  between  the  electrode 
potentials  and  overpotentials  can  be  found  in  Appendix  B. 

After  separating  overpotentials  of  rj3pB,  r\2 PB  and  r/s  from  respec¬ 
tive  electrode  potentials,  Eqs.  (4)  and  (6)-(8)  can  be  rewritten  into 
Butler-Volmer  type  expressions  [28,39]. 

rS2  =  rs 2,0  ( exp(-as//?s)  -  ,ad  exp[(l  -  as^s]}  (14) 

[^  L0  ,ad,eq  J 


r S4  =  ts4,o  <  ^ — —  cxp(  —0i3pBfv)3pB )  -  exp[(l  -  a3pB)frj3pB]  > 


(15) 


f  Or  ^V.MIEC 

rB3  =  rB3,o  <  r - ^~r - exp(asfrjs)  -  exp[-(l  -  ots)fris] 

L0~  eqCV,MIEC,eq 


(16) 


rB4  =  1~B4,0 


exp{-2oi2PBfri2PB) 


Cy.MIEC 

CyMIEC.eq 


exp[2(l  -  ^PBlf^PB] 


(17) 


By  setting  the  forward  and  backward  reaction  rates 
equal  to  each  other  at  OCV,  the  exchange  reaction  rates 
rs2,o,  rS4,o,  rB3,o  and  rB4,o  can  be  derived  by  replacing  E0cv 
with  constants  as: 


rs2,o  =  k< 


S2 


ks2fcsiTp(02)i/2 

KS2KS1 


(l-«s) 


(C0-.. 


rs4,0  =  1<S4 


kS3  (fes4)1_“3PB 

k-3  (ks-4cv,Yszr)(-“3PBi 


(C0  ,,eq) 


l-a3pB 


(18) 

(19) 


rB3,0  =  Cfes3 J*** CfeB3 )"*  “S(/TS(C0-  eqCV.MIEC.eq)1  (20) 

ad’  M 

tB4,0  =  (^B4Cv,YSz)1_a2PB(^B4Cv,MIEC,eq)a2PB  (21) 


3.  Formulation  of  the  1-D  numerical  model 


requires  that  specific  RLSs  be  screened  out  for  simplifying  the  cal¬ 
culation  requirements  and  relating  the  simulation  to  experimental 
findings.  As  for  a  LSM-type  cathode  with  low  intrinsic  ionic  conduc¬ 
tivity,  oxygen  exchange  reactions  involving  an  oxygen  vacancy  (B3 
and  B4)  can  be  regarded  as  rate-limiting,  van  Heuveln  and  Jiang’s 
works  [10,17]  also  identified  surface  electron  donation  (S2)  and 
oxygen  diffusion  (S3)  as  the  most  possible  RLSs  for  the  3PB  kinet¬ 
ics  on  LSM  cathode.  For  this  reason,  steps  S2,  S3,  B3  and  B4  are 
set  as  RLSs  with  finite  net  rates,  and  the  steps  SI  and  S4  remain  at 
quasi-equilibrium  with  net  rates  ignored.  This  scenario  highlights 
contributions  of  electrochemical  reactions  to  the  transition  of  bi¬ 
pathway  cathode  kinetics  in  this  work,  which  is  the  main  difference 
from  previous  models  [33-35].  The  coupled  mass-transport  equa¬ 
tions  can  be  built  for  the  two  dependent  variable  species  using  Fields 
2nd  law  for  diffusion: 


9C0-,a 

dt 


n  /92C0-,ad\  ,  ,  ,  n  /^Qr.adN 

Ds,chem  ^  ^2  J  +  ^S2  ^  “  Ds,chem  ^  ^  J 

+  rS2,o  ( exp(-asfrjs)  -  C°  ’ad  exp[(l  -  o?sfe]} 

L  L0-,ad,eq  J 


C0  CV.MIEC 


9Cv.iv 


dt 


=  £>b, 


Cq-^  j  eq  Cv.MIEC.eq 
f  92Cy,MIEC 


exp(o'^s)  -  exp[-(l  -  as)frjs] 


} 


(22) 


,hem 

AS 


AS  (  d2Cy,MiEC  N 

^7rB3-Db,chem  ^  ax2  ^ 


AV 


^B3,0 


L  L0-,ad,< 


addv.MIEC 


,eq  CV.MIEC.eq 


exp(as/^s)-exp[-(l  - 


as)/>7s]j  (23) 


where  Dschem  and  Db  chem  are  the  surface  and  bulk  chemical  diffu- 
sivities  of  oxygen  and  oxygen  vacancy,  respectively;  and  AS/ A  V  is 
the  volume-specific  surface  area. 

At  the  MIEC/EC  boundary  (x  =  lc  in  Fig.  2),  a  zero-flux  boundary 
condition  is  considered  for  surface  oxygen  ion  as  in  Eq.  (24),  while 
for  bulk  path  a  blocking  boundary  condition  is  effectively  estab¬ 
lished  with  Eq.  (25)  assuming  oxygen  vacancy  concentration  at  this 
boundary  should  be  same  as  that  at  equilibrium,  in  accordance  to 
the  physical  assumptions  in  Fig.  2: 

AS  r>  f^O~,ad  . 

^yf^s.chem  (  )  -  rS2  rB3 


dx 


x=\c 


=  rs 2,0  <  exp(-asfr]s) - 


C°a"d’eCl 


exp[(l  -  as)/r/s 


4*  Db.chem 


9CV, 


,MIEC 


dx 


(24) 


x=lc 


Cv,MIEC,x=lc  =  ^v.MIEC.eq  (25) 

At  the  MIEC/electrolyte  interface  (x  =  0),  with  rS4  ignored  at  equi¬ 
librium  it  can  be  derived  that 

C0-  =  C0-  eq  exp(/773PB)  (26) 

ad  ad’  M 


t^b.chem 


Vx=0  =  _rB4  =  ^4,0 


CV.MIEC 

Cy.MIEC.eq 


exp[2(l  -  Q(2Pb)/^2Pb] 


-  exp(-2o'2PB/^2PB) 


(27) 


3.1.  Governing  equations  and  boundary  conditions 


The  model  considers  C0-  and  Cv.miec  as  two  independent  vari- 

ad 

ables.  Concentrations  of  other  species  are  taken  as  constant  or 
simple  functions  of  the  two  primary  variables.  This  condition 


Flence,  we  have  a  Dirichlet-type  boundary  condition  for  the  sur¬ 
face  oxygen  ions  arising  from  a  charge-transfer  reaction  assumed 
fast-enough  at  3PB,  and  a  Neumann-type  boundary  condition  for 
oxygen  vacancies  of  the  MIEC  bulk  assuming  a  relatively  slow  2PB 
exchange  reaction. 


M.  Gong  et  al.  /  Journal  of  Power  Sources  201  (2012)  204-218 


209 


Table  2 

Values  of  type  I  and  II  parameters  used  in  the  simulation  cases. 


Parameter 

Case  1 

Case  2 

Case  3 

Description 

f^b.chem 

1  x  10-6  cm2  s_1 

Same  (as  Case  1) 

Same 

Surface  oxygen  diffusivity 

Hs.chem 

1  x  10-6  cm2  s  1 

Same 

Same 

Bulk  vacancy  diffusivity 

Cv.MIEC.eq 

1  x  10-8  mol  cm-3 

1  x  10-8  mol  cm-3 

1  x  10-7  mol  cm-3 

Equilibrium  vacancy  cone,  in  MIEC 

C°ad’eq 

1  x  10-11  molarr2 

1  x  10_1°  mol  cm-2 

1  x  10_1°  mol  cm-2 

Equilibrium  surface  oxygen  cone. 

fS2,0 

1  x  10-7  mol  cm-2  s-1 

Same 

Same 

Exchange  rate  constant  for  S2 

053,0 

1  x  10-7  mol  cm-2  s-1 

Same 

Same 

Exchange  rate  constant  for  B3 

054,0 

1  x  10-6  mol  cm-2  s_1 

Same 

Same 

Exchange  rate  constant  for  B4 

AS(AV)-1 

1  x  105  cm-1 

Same 

Same 

Volume-specific  surface  area 

II  Oi 3PB,  C^PB,  OLs 

0.5 

Same 

Same 

Symmetry  factor 

r 

10-9  mol  cm-2 

5  x  10-9  mol  cm-2 

5  x  10-9  mol  cm-2 

Surface  adsorption  site  density 

Oeq 

0.01 

0.02 

0.02 

Equilibrium  oxygen  coverage 

T 

1073  K 

Same 

Same 

Temperature 

p02 

0.21  atm 

Same 

Same 

Oxygen  partial  pressure 

3.2.  Numerical  method 

Instead  of  direct  evaluation  of  the  steady-state  solutions  as 
in  previous  bi-pathway  models,  a  finite  control-volume  analy¬ 
sis  was  adopted  together  with  time-discretization  to  obtain  the 
transient  solution  of  the  variables,  thus  allowing  elucidation  of 
non-stationary  electrode  kinetics.  The  modeled  distance  (10|jim) 
from  the  MIEC/electrolyte  interface  (x  =  0)  to  MIEC/EC  bound¬ 
ary  (x  =  lc)  is  divided  into  40  flux  nodes  with  equal  spacing. 
The  locations  for  node  properties  (e.g.  specie  concentration)  are 
offset  by  1/2  grid  spacing  away  from  the  flux  nodes  following 
Patankar’s  method  with  staggered  node  matrix  [41].  The  com¬ 
putational  code  is  implemented  with  Visual  C++  6.0  software. 
Overpotentials  of  +0.1  V  to  -0.45  V  are  given  as  the  external  input 
with  50-100  mV  intervals  between  data  points.  At  t= 0  computa¬ 
tion  starts  with  initialization  of  surface  and  bulk  species  profiles 
using  input  boundary  conditions,  and  attains  steady-state  after 
a  calculation  period  of  240  nominal  time-steps,  which  actually 
spans  from  a  few  hours  to  tens  of  hours  according  to  input  val¬ 
ues  and  material  parameters.  Calculation  of  the  model  requires 
a  large  number  of  physical  parameters,  which  are  roughly  classi¬ 
fied  into  type  I  and  II.  Type  I  are  measurable  material  parameters 
directly  used  in  governing  equations  and  boundary  conditions  for 


simulation.  Type  II  includes  material  parameters  like  the  defined 
forward/backward  rate  constants  not  directly  used  in  numerical 
solution,  and  other  thermal  and  electrochemical  condition  con¬ 
stants.  The  parameter  values  for  each  case  are  listed  in  Table  2. 
Type  I  parameter  values  have  been  given  in  previous  bi-pathway 
models  except  for  C0-  [33,35],  but  were  more  suited  for  sim- 

ulating  MIECs  with  good  ionic  conductivity.  Re-estimation  of  the 
values  is  thus  made  from  LSM  properties  established  by  iso¬ 
topic  oxygen  exchange  (IOE)  and  electrical  conductivity  relaxation 
(ECR)  measurements  [42-46].  The  parameter  values  in  Table  2 
are  believed  to  be  within  order  of  magnitude  accuracy  at  best 
due  to  the  lack  of  experimental  data  for  the  simulated  condi¬ 
tions.  Details  on  parameter  selection  are  separately  discussed  in 
Appendix  A. 


4.  Results  and  discussion 

The  simulation  was  performed  with  governing  Eqs.  (22)  and 
(23),  boundary  conditions  (24)-(27)  and  Table  2  parameters.  The 
transient  profiles  for  electrode  species  are  given  in  Section  4.1  to 
monitor  the  convergence  and  efficiency  of  the  model.  Steady-state 


<2 

cr 

c 

TT 

< 

£D 

O 

D 

3 

O 

*< 

o 

o 

3 

o 

ft) 

=3 


D 

r-+ 

o' 

3 

o 


3 

o 

• 

o 

3 


Fig.  3.  Transient  diffusion  profiles  of  surface  oxygen  ion  and  bulk  vacancy  calculated  with  CV.eq  =  1  x  10-7  mol  cm-3  and  C0-  eq  =  1x10  10  mol  cm  2  at  -0.2  V  overpotential 
(open  symbol:  C0- ;  filled  symbol:  CV;  value  of  actual  time-step:  3.08  x  10-9  s;  series  1  =2.17  x  10-3  s). 
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Fig.  4.  Transient  profiles  of  the  bulk  vacancy  flux  calculated  with 
Cv,eq  =  l  x  10_7molcm“3  and  C0-eq  =  1  x  10_10mol  cm-2  at  -0.2  V  overpotential 
(series  1  =2.17  x  10-3  s). 

results  are  discussed  in  Sections  4.2  and  4.3  to  investigate  kinetic 
transition  of  the  MIEC  cathode  from  case  to  case. 

4. t.  Transient  distribution  of  the  active  electrode  species 

Figs.  3  and  4  show  the  transient  profiles  for  concentration  and 
diffusion  flux  of  surface  oxygen  ion  (C0- )  and  bulk  oxygen  vacancy 
(Cy)  under  moderate  overpotential  (-0.2  V).  Values  of  C0-  and 

Cv.MiEc.eq  are  chosen  as  those  reported  for  LSM  material.  A  cathodic 
overpotential  of  -0.2  V  is  applied  to  represent  typical  SOFC  oper¬ 
ation  condition.  Transient-state  calculation  increases  this  model’s 
applicability  compared  to  previous  steady-state  analysis  [33-35], 
as  the  data  can  be  used  to  predict  AC  impedance  evolution  during 
electrode  activation  and  also  assess  the  simulation  efficiency.  The 
actual  time-step  used  for  simulation  is  about  3  x  10-9  s,  which  is 
chosen  to  meet  the  stability  and  plausibility  criteria  from  Patankar 


Fig.  5.  Logarithmic  current  versus  overpotential  profiles  for  bi-pathway  constituted 
overall  electrode  kinetics  with  Case  1  parameters  as  C0-  eq  =  1  x  10-11  mol  cm-2 
and  Cv,eq  =  1  x  10-8  mol  cm-3. 

[41  ].  This  however,  makes  convergence  slow.  To  observe  how  the 
simulation  reaches  steady-state,  transient  profiles  are  plotted  out 
at  certain  stages  of  the  actual  simulation  process,  which  are  denoted 
as  series  1-240  to  represent  the  increment. 

Fig.  3  displays  surface  oxygen  concentration  distributes  in  a  way 
reverse  to  the  bulk  oxygen  vacancy,  and  their  profiles  reflect  the 
set-up  of  reaction  rates  as  boundary  conditions  at  the  two  inter¬ 
faces.  The  species  concentrations  in  Fig.  3  seem  to  stabilize  fast,  but 
Fig.  4  analysis  shows  the  adjustment  of  the  bulk  vacancy  flux  near 
the  MIEC/EC  interface  lasts  until  simulation  ends.  This  late  conver¬ 
gence  at  diffusion  front  of  oxygen  vacancy  spans  up  to  tens  of  hours. 
The  electrode  is  further  activated  at  regions  adjacent  to  x  =  0  (by 
increase  of  vacancy  flux),  and  less  activated  beyond  2-3  p,m  away 
from  the  MIEC/electrolyte  interface  where  2PB  current  decays  to 


Fig.  6.  Evolution  of  the  local  reaction  current  (i'ioc)  in  the  total  3PB  current  (it0 1)  as  function  of  overpotential  in  Case  1  with  C0-  eq 
Cv.eq  =  1  x  10-8  mol  cm-3. 


=  1x10  11  mol  cm  2 


and 
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Fig.  7.  Logarithmic  current  versus  overpotential  profiles  with  parameters  in  Case  2 
as  C0-  eq  =  1  x  10-10  mol  cm-2  and  Cv,eq  =  1  x  10“8  molcirr3. 


mA  range.  Analysis  on  steady-state  species  distribution  and  active 
reaction  zone  as  functions  of  overpotential  and  material  properties 
is  done  in  Section  4.3. 


4.2.  Kinetic  transition  of  the  dominant  pathway 

Relative  dominance  of  the  3PB  and  2PB  pathways  at  different 
overpotentials  is  investigated  to  understand  how  the  kinetic  transi¬ 
tion  is  affected  by  material  parameters  and  its  consequent  impacts 
on  electrode  performances.  Figs.  5,  7  and  8  are  the  logarithmic 
current-overpotential  (/-V)  curves  with  separate  kinetic  pathway 
contributions  for  Case  1,  2  and  3,  respectively.  The  simulation  is 
performed  between  -0.4  V  to  +0.1  V  with  50-100  mV  intervals,  and 
the  equilibrium  concentrations  of  surface  oxygen  ion  (C0-  eq)  and 
bulk  oxygen  vacancy  (Cv,eq)  are  varied  among  the  three  cases  for  a 
parametric  study  on  electrode  kinetics: 


(1)  C0-  eq  =  1.0  x  10  11  mol  cm  2 
for  Case  1 

(2)  C0-  eq  =  1.0  x  10_10mol  cm”2 
for  Case  2 

(3)  C0-  eq  =  1.0  x  10_10mol  cm”2 
for  Case  3 


&  Cv,eq  =  1.0  x  10-8  mol  cm-3 
&  Cv,e q  =  1.0x  10_8molcnrr3 
&  Cv,eq  =  1.0  x  10-7  mol  cm-3 


Lines  are  drawn  in  the  figures  to  display  the  Tafel  estimation 
of  the  exchange  current  densities  for  the  two  paths.  The  exchange 
current  densities,  kinetic  transition  voltage  and  overall  current  den¬ 
sities  at  -0.3  V  identified  for  each  case  are  listed  in  Table  3.  The  3PB 


Fig.  8.  Logarithmic  current  versus  overpotential  profiles  with  parameters  in  Case  3 
as  C0-  eq  =  1  x  lO-10  mol  cm-2  and  Cv,eq  =  1  x  10-7  mol  cm-3. 


and  2PB  current  densities  are  given  by  the  fluxes  of  surface  oxygen 
ions  and  bulk  oxygen  vacancies  at  x  =  0  respectively  as: 


■A£i 

AV 


*3PB  —  2F  a  ^Ds.chem 


dC0- 


,ad 


dx 


(28) 


*2PB  —  ^FDb.chem 


f  9Cv,miec  \ 

V  9*  ) 


x=0 


(29) 


4.2 A.  Overpotential  effects  on  pathway  transition 

Simulation  results  from  the  I-V  profiles  and  corresponding 
exchange  current  densities  are  used  to  analyze  overpotential  effects 
on  the  kinetic  pathway  transition  for  MIEC  cathode  with  Case  1 
parameters.  With  low  values  of  both  C0-  eq  and  Cv,eq  in  Case  1, 
Fig.  5  demonstrated  that  the  3PB  pathway  dominates  at  relatively 
low  cathodic  and  anodic  overpotentials,  while  the  2PB  pathway 
current  become  increasingly  important  at  higher  cathodic  overpo¬ 
tentials.  A  change  of  dominance  on  electrode  kinetics  from  3PB  path 
to  2PB  path  is  observed  at  around  -0.33  V,  resulting  in  a  deviation 
of  the  total  current  density  from  the  3PB  current  profile  and  an 
improved  electrode  performance.  Since  this  transition  of  electrode 
kinetics  at  negative  overpotential  is  the  main  research  interest,  no 
further  simulation  in  the  anodic  polarization  region  (>0.1  V)  was 
made. 

The  logarithmic  relationship  between  current  density  and  over¬ 
potential  in  Fig.  5  enables  a  Tafel  analysis.  Currents  from  both 
pathways  enter  into  the  Tafel  zone  at  overpotential  around  -0.1  V, 
close  to  value  given  by  high-field  approximation  of  the  B-V  equa¬ 
tion  at  800 °C  (-fy»RT/F^  0.09  V)  [10,47].  Flowever,  the  dashed 
tangent  lines  of  the  profiles  yield  rather  different  exchange  cur¬ 
rent  densities  (z’o).  For  the  3PB  pathway,  the  z'o,3pb  intercept  is  about 
1  o-2  05  A  cm-2,  matching  well  with  the  zS2,o  calculated  from  rS2,o 
for  surface  electron  donation.  Tafel  analysis  in  Fig.  5  finds  a  rather 


Table  3 

Summary  of  V-I  curve  simulaiton  and  tafel  analysis  results  for  Cases  1-3. 


Case 

Cv.eq 

(mol  cm-3) 

C0  ,eq 

(mol  cm-2) 

Transition 
voltage  (V) 

Set-up  of 

lO,3PB 

(mAcirr2) 

Extracted 

lO,3PB 

(mA  cm-2) 

Set-up  of 

io,2PB 

(mAcirr2) 

Extracted 

lO,2PB 

(mAcirr2) 

Extracted 

lO.TOT 

(mAcirr2) 

i'tot  at 
-0.3  V 
(mAcirr2) 

Control  of 
electrode 
activation 

1 

1  x  10“8 

1  xlO”11 

-0.33 

10 

10 

193 

7x  10“5 

7 

163 

3PB & 2PB 

reactions 

2 

1  x  lO”8 

1  xlO-10 

-0.45 

10 

49 

193 

00 

X 

o 

45 

598 

Diffusion 

3 

1  x  10“7 

1  xlO”10 

-0.26 

10 

45 

193 

7x  10“4 

45 

719 

Diffusion 

212 


M.  Gong  et  al.  /  Journal  of  Power  Sources  201  (2012)  204-218 


small  2PB  exchange  current  density  z0,2pb  of  10-417  A  cm-2,  com¬ 
paring  to  an  input  rB 4,0  equivalent  to  about  1 0-2  A  cm-2.  The  overall 
exchange  current  density  is  close  to  z0,3pB,  around  6-7mAcm-2. 
The  results  indicate  that  the  electrode  activation  in  Case  1  is 
mainly  controlled  by  local  reaction  S2  to  form  intermediate  O-. 
The  large  deviation  of  z0,2pb  from  the  input  value  implies  2PB  kinet¬ 
ics  is  less  favored  by  factors  other  than  electrochemical  force.  At 
lower  overpotential,  electrochemical  terms  in  Eqs.  (14)-(17)  are 
less  influential,  thus  O-  from  step  S2  and  oxygen  vacancy  from  2PB 
exchange  step  B4  could  chemically  favor  a  positive  direction  of  step 
B3  reaction,  slowing  down  the  increase  of  vacancy  concentration 
and  flux.  Additionally  the  volume-specific  surface  area  ( AS/AV)  as 
a  geometric  factor  in  Eq.  (28)  also  introduces  higher  3PB  current. 
Hence  a  larger  driving  force  for  vacancy  exchange  is  needed  from 
more  negative  overpotentials,  so  that  the  2PB  reaction  can  out¬ 
weigh  the  combined  geometry  and  oxygen  incorporation  effects  to 
dominate  over  the  3PB  pathway  kinetics.  As  a  result  of  the  kinetic 
transition  to  slower  2PB  kinetics,  the  total  exchange  current  den¬ 
sity  exhibits  a  rather  lower  value  than  the  3PB  exchange  current 
density  in  this  case. 

4.2.2.  The  mechanism  of  pathway  transition 

To  further  understand  the  interaction  between  3PB  and  2PB 
pathways  during  electrode  activation,  the  composition  of  3PB  cur¬ 
rent  is  analyzed  in  terms  of  diffusion  and  local  reaction  currents  to 
unveil  the  mechanistic  process  potentially  involved  in  the  pathway 
transition.  In  Figs.  6  and  9  the  change  of  local  current  density  zjoc  as 
a  function  of  overpotential  is  examined  for  each  case.  The  total  3PB 
current  density  zt0 1  can  be  viewed  as  the  sum  of  diffusion  current 
density  /diff  from  the  accumulative  flux  of  O-  just  outside  TPB  and 
local  current  density  zjoc  from  surface  reactions  rs 2  and  rB3  at  x  =  0, 
and  then  idiff  can  be  expressed  as 

Miff  =  *3PB,tot  -  *loc  (30) 

zjoc  can  be  derived  from  the  two  RLSs  as  [10,48] 


general  the  analysis  on  3PB  current  composition  suggests  limita¬ 
tion  of  surface  oxygen  diffusion  as  a  potential  mechanistic  process 
by  which  kinetic  pathway  transition  occurs.  This  transition  mech¬ 
anism,  however,  has  seldom  been  reported,  van  Heuveln  found  no 
diffusion-limiting  current  for  cathode  kinetics  governed  by  surface 
pathway  within  -0.3  V  overpotential  [10],  while  at  more  negative 
overpotential  Siebert  found  impedances  typical  of  mass-transport 
limiting  process  during  cathode  kinetics  transition  to  bulk  path 
dominance  [16].  In  addition,  Coffey’s  simulation  also  showed  a 
diffusion-limit  for  3PB  current  after  the  pathway  transition,  though 
no  specific  process  was  identified  [35].  Comparing  with  the  last 
results,  it  is  clarified  in  this  model  that  limiting  of  surface  oxygen 
diffusion  is  the  physical  process  behind  surface-to-bulk  pathway 
transition. 

4.2.3.  Effect  of  material  property  on  pathway  transition 

Figs.  7  and  8  show  the  effects  from  change  of  material  prop¬ 
erty  parameters  C0-  eq  and  Cv>eq  on  the  kinetic  pathway  transition. 
In  Case  2,  where  a  higher  value  of  C0-  eq  is  used  while  the  other 
parameters  remain  the  same  as  Case  1,  3PB-to-2PB  kinetic  transi¬ 
tion  is  postponed  more  negatively  to  overpotential  around  -0.45  V 
and  current  above  1 A  cm-2.  In  contrast,  Fig.  8  shows  that  the  transi¬ 
tion  of  the  dominant  pathway  occurs  at  more  positive  overpotential 
of  around  -0.26  V  due  to  the  increase  of  Cv,eq  by  10  times  in  Case 
3  relative  to  Case  2.  The  reasons  for  such  kinetic  discrepancy  for 
different  cases  are  explored  with  Tafel  analysis  and  3PB  current 
division. 

Tafel  analysis  in  Fig.  7  yields  z0,3pB  of  10-1,42  Acnrr* 1 2  and  z0,2pb 
of  10-409  A  cm-2  for  Case  2.  The  overall  exchange  current  density 
is  almost  identical  to  the  3PB  one,  and  is  not  reduced  by  sluggish 
2PB  activation  as  in  Case  1.  z0,3pB  obtained  from  Fig.  7  is  consider¬ 
ably  higher  than  calculated  with  input  S2  and  B3  rate  constants. 
This  result  may  suggest  surface  reactions  of  O-  no  longer  dominate 
3PB  path  kinetics.  As  indicated  by  Fig.  9,  the  local  current  density 
estimated  with  Eqs.  (30)  and  (31)  accounts  for  less  than  20%  ofthe 


Moc  —  *anodic,loc  *cathodic,loc 


iB35o(cv.MiEc,x=o/Cv,MiEc,eq)exp[(2ci's)/^s]  •  iS2,o{exp[(2  -  as)frjs]  -  exp(-ajT]s)} 

{iS2,0  exp[(2  -  as)fr]s]  +  iB3,o(Cv,MiEc,x=o/Cv,MiEc,eq)exp[(l  +  as)fr]s]}  •  feo  exp(-asfr)s )  +  iB3,o  exp[-(l  -  as)fr]s]} 

ij2  0  ■  lB3,0  exp[2(l  -  0's)/^s]{(Cv,MIEC)x=o/CV,MIEC,eq)exp[(l  +  as)fijs]  -  exp[-(l  -  as)fijs]} 

{is2,o  exp[(2  -as)fr]s]  +  iB3,o(CvjviiEc,x=o/CvjviiEc,eq)exp[(l  +as)fr]s]}  ■  {is2,o  exp(-Q's/^s)  +  i'b3,o  exp[-(l  -as)fils]} 


(31) 


Case  1  study  in  Fig.  6  shows  that  at  -0.1  to  -0.2  V  most  of  the 
3PB  current  density  comes  from  the  local  reaction  current,  con¬ 
firming  that  electrode  activation  is  mainly  controlled  by  surface 
reactions  of  rS2  and  rB3.  The  ratio  between  zjoc  and  idiff  decreases 
as  cathodic  overpotential  increases,  and  reaches  a  minimum  just 
before  the  3PB-to-2PB  transition  at  -0.3  to  -0.4  V.  The  increas¬ 
ing  contribution  from  diffusion  current  at  -0.1  to  -0.3  V  indicates 
that  electrode  surface  beyond  3PB  becomes  more  active  for  oxygen 
reduction  during  the  polarization,  and  the  effective  reaction  zone 
is  expected  to  extend  accordingly  towards  the  MIEC/EC  boundary. 

Figs.  5  and  6  suggest  by  what  mechanism  the  electrode  kinet¬ 
ics  switch  from  one  dominant  pathway  to  the  other.  The  total  3PB 
current  in  Case  1  encounters  a  diffusion  limit  associated  with  sur¬ 
face  oxygen  flux  above  -0.3  V,  suggested  by  the  increase  of  zjoc/itot- 
The  coincidence  of  the  diffusion  limit  for  3PB  current  with  kinetic 
transition  to  2PB  dominance  implies  then  that  the  bulk  vacancy 
concentration  increases  faster  than  the  production  of  O-  via  rs 2 
outside  of  the  3PB.  Consequently,  additional  surface  oxygen  diffu¬ 
sion  flux  towards  3PB  is  prevented  by  local  consumption  via  rB3, 
e.g.  a  surface  diffusion  limit  occurs.  Meanwhile,  local  reaction  cur¬ 
rent  (triangle  in  Fig.  6)  at  3PB  appears  to  be  less  limited,  and  with 
a  linear  overpotential-dependence.  This  is  because  the  concentra¬ 
tion  of  O-  at  x  =  0  is  directly  defined  by  the  boundary  condition, 
and  the  driving  forces  for  reactions  S2  and  B3  are  almost  certain.  In 


total  3PB  current  for  Case  2,  implying  surface  oxygen  diffusion  is 
now  the  dominant  3PB  step.  Due  to  the  change  of  electrode  activa¬ 
tion  to  surface  diffusion  control  at  higher  C0-  eq,  both  3PB  kinetics 
and  the  overall  electrode  performance  improve  compared  to  Case 

1.  Increase  of  C0-  eq  in  Case  2  would  alter  the  absolute  concen¬ 
tration  of  O-,  though  C0-/C0-  eq  is  unchanged  at  3PB  for  certain 
overpotential.  Consequently  this  introduces  a  rather  higher  level 
of  surface  oxygen  diffusion  flux  at  3PB  in  Fig.  9,  and  gives  rise  to  a 
larger  overall  3PB  current.  To  overcome  the  increased  suppression 
effects  from  3PB  path  at  high  C0-  eq,  the  oxygen  vacancy  concen¬ 
tration  at  2PB  has  to  rise  correspondingly  higher,  and  hence  the 
overpotential  for  pathway  transition  grows  more  negative  in  Case 

2.  As  for  Case  3,  the  higher  Cv,eq  used  in  simulation  brings  about  a 
reduction  of  the  total  3PB  current  and  concomitant  rise  of  zjoc/*tot 
ratio  in  Fig.  9.  Basically,  this  change  reflects  increasing  contribution 
from  2PB  kinetics  with  the  higher  value  of  Cv,eq,  as  both  concentra¬ 
tion  and  flux  of  bulk  vacancies  would  increase  accordingly.  The  2PB 
exchange  current  for  Case  3  also  increases  in  Fig.  8  at  higher  Cv,eq, 
resulting  in  a  slightly  decreased  overall  exchange  current  densities 
comparing  to  z0>3pB. 

Notably,  this  observation  reveals  that  changes  in  surface  and 
bulk  parameters  affect  electrode  kinetics  in  different  manners.  An 
increase  of  oxygen  vacancy  equilibrium  concentration  promotes 
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Fig.  9.  Comparison  of  the  relationship  between  local  reaction  currents  and  total  3PB  currents  in  Case  2  and  Case  3  (Case  2:  C0- eq  =  1  x  10  10  mol  cm  2  and 
Cv>eq  =  1  x  1 0“8  mol  cm"3 ;  Case  3 :  C0-  eq  =  1  x  1 0“10  mol  cm-2  and  Cv,eq  =  1  x  1 0~7  mol  cm"3 ). 


an  early  transition  to  bulk  pathway  and  actually  limits  the  oxy¬ 
gen  surface  diffusion,  while  large  equilibrium  coverage  of  surface 
oxygen  ion  improves  3PB  kinetics  by  enhancing  oxygen  transport 
on  the  surface.  Despite  the  differences  of  electrode  kinetics  and 
performance  between  Cases  1-3,  Figs.  6  and  9  did  show  that  the 
response  of  i\oclkot  ratio  to  overpotential  follows  a  similar  trend  for 
all  the  cases  studied,  confirming  the  limitation  of  surface  diffusion 
of  oxygen  ion  as  the  physical  mechanism  for  3PB-to-2PB  pathway 
transition. 

4.2.4.  Association  to  electrode  performance  and  experimental 
finding 

Combined  simulation  results  from  Cases  1-3  indicate  that 
the  transition  of  the  dominant  pathway  greatly  depends  upon 
material  properties,  with  earlier  transition  in  terms  of  cathodic 
overpotential  promoted  by  decreasing  equilibrium  surface  oxy¬ 
gen  concentration  C0-  eq  and  increasing  equilibrium  bulk  oxygen 
vacancy  concentration  Cv,eq-  Change  of  C0-  eq  and  Cv,eq  also 
determines  the  RLS  for  electrode  kinetics,  and  accordingly  the 
performance  in  Table  3.  From  the  practical  viewpoint,  lower  over¬ 
potential  transition  of  both  2PB  and  3PB  paths  would  be  preferred 
to  increase  total  current  density.  Flowever,  from  Table  3  it  can 
be  learned  that  earlier  3PB-to-2PB  transition  does  not  necessarily 
lead  to  higher  electrochemical  performance.  Case  3  has  the  highest 
current  density  at  -0.3  V,  due  to  the  increased  2PB  kinetic  contri¬ 
bution  at  a  relatively  smaller  overpotential.  But  the  dominance  of 
3PB  kinetics  in  Case  2  leads  to  greater  performance  improvement 
compared  to  Case  1,  even  with  a  late  transition  voltage.  Generally, 
the  value  of  C0-  eq  can  be  viewed  to  associate  with  electrocatalytic 
activity  of  the  cathode  material,  while  Cv,eq  relates  to  ionic  conduc¬ 
tivity  for  the  MIEC.  For  practical  development  of  high  performance 
SOFC  cathode,  usually  poor  MIEC  like  LSM  is  either  directly  mixed 
with  a  good  ionic  conductor  like  YSZ  or  impregnated  with  ionic 
conductive  and  catalytic  active  materials,  such  as  doped  ceria.  The 
simulation  scenario  in  this  study  imply  that  improving  surface  cat¬ 
alytic  activity  for  a  poor  MIEC  with  certain  ionic  conductivity  can 
be  more  efficient  in  terms  of  electrode  performance  promotion  at 
a  defined  electrode  configuration  (fixed  3PB  and  2PB  contact  area). 


In  comparison  to  the  simulation  results,  it  is  well  known  for 
experimental  and  theoretical  studies  that  SOFC  cathode  kinetics 
can  be  improved  at  cathodic  overpotentials  promoting  oxygen 
reduction  through  bulk  pathway.  However,  the  reported  overpo¬ 
tential  when  2PB  contribution  becomes  apparent  tends  to  vary 
from  case  to  case.  Siebiert  noticed  considerable  current  increase  for 
LSM  cathodes  around  overpotential  of  -0.4  V  [16].  Kim  reported 
the  resistance  of  LSM-YSZ  electrode  was  greatly  reduced  after 
current  treatment  at  -0.5  V  [12].  At  higher  temperature  (950 °C), 
van  Heuveln  found  significant  relaxation  of  LSM  cathode  resis¬ 
tance  shortly  after  DC  polarization  at  -0.3  V  [10].  More  recently, 
studies  from  Fleig’s  group  and  Shao-Horn’s  group  on  dense  LSM 
microelectrodes  at  700-800 °C  both  indicated  that  the  bulk  path 
contribution  increased  dramatically  at  DC  bias  of -300  mV  [22,24]. 
In  this  work,  the  simulated  overpotential  for  the  pathway  transition 
is  around  -0.26  V  to  -0.45  V,  lying  within  the  reported  data  range. 
The  result  discrepancy  rather  reflects  that  the  kinetic  transition 
between  surface  and  bulk  path  for  oxygen  reduction  is  sensitive 
to  electrode  operation  conditions  and  material  properties,  which 
closely  depend  on  the  manufacturing,  composition,  geometry  and 
structure  of  the  electrode. 


4.3.  Dependence  of  active  reaction  zone  and  surface  reaction 
rates  on  overpotential  and  material  parameters 

In  the  physical  scenario  of  oxygen  reduction  (Fig.  2),  a  “critical” 
boundary  is  assumed  to  exist  at  x  =  lc  beyond  which  surface  and 
bulk  exchange  reactions  would  cease  to  modify  the  local  vacancy 
concentration  from  Cv,eq.  The  active  reaction  zone  thickness  repre¬ 
sents  the  electrode  area  actually  participating  in  oxygen  reduction, 
and  is  important  to  cathode  performance  improvement  through 
microstructure  design.  A  detailed  analysis  is  thus  performed  in  Sec¬ 
tions  4.3.1  and  4.3.2  to  derive  the  thickness  of  active  reaction  zone 
from  species  concentration  profiles  and  identify  its  dependence  on 
overpotential  and  material  property  from  reaction  rates  rs 2  and  rB 3. 
The  results  are  used  to  relate  simulation  findings  to  the  practical 
routes  for  cathode  performance  improvement. 
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Electrode  thickness(nm) 

Fig.  10.  Distribution  profiles  of  oxygen  vacancy  concentrations  under 
different  overpotential  in  Case  1  as  C0-  eq  =  1  x  10-11  mol  cm-2  and 
Cv.eq  =  1  x  10-8  mol  cm-3. 

43  A.  Effects  of  overpotential  and  material  property  on  active 
reaction  zone 

Figs.  10  and  11  show  the  overpotential-driven  profiles  of  the 
normalized  oxygen  vacancy  concentration  along  the  electrode  in 
Case  1  and  Case  2,  respectively.  The  electrode  is  electrochemi- 
cally  activated  near  2PB/3PB  interfaces  as  indicated  by  increased 
local  vacancy  concentrations.  Applying  more  negative  overpoten¬ 
tial  increases  the  penetration  depth  of  oxygen  vacancy  deviating 
from  thermal  equilibrium  concentration.  In  the  magnified  views  of 
Figs.  10  and  11,  the  active  reaction  zone  thickness  can  be  seen  to 
extend  from  ~1  |jim  in  Case  1  to  ~4  p,m  in  Case  2  at  higher  value 
of  C0-  .  Accordingly,  local  vacancy  concentrations  outside  of  2PB 

in  Case  2  also  increase  from  Case  1  results.  This  suggested  change 
of  surface  material  properties  not  only  affects  3PB  pathway  trans¬ 
port  but  also  the  bulk  pathway  kinetics  via  a  vacancy  concentration 
adjustment.  This  phenomenon  will  be  addressed  in  later  discussion 
on  methodology  for  cathode  performance  improvement.  The  rea¬ 
son  for  extension  of  the  active  reaction  zone  is  discussed  with  Fig.  1 2 


Fig.  11.  Distribution  profiles  of  oxygen  vacancy  concentrations  under 
different  overpotential  in  Case  2  as  C0- eq  =  1  x  10-10  mol  cm ~2  and 
Cv.eq  =  1  x  10-8  mol  cm-3. 


results.  An  increase  of  C0-  eq  by  10  times  in  Case  2  relative  to  Case 
1  reduces  the  local  ratio  of  C0-/C0-  eq  in  Fig.  12  while  increasing 
the  surface/bulk  exchange  reaction  rates  rB 3.  It  is  noteworthy  that 
such  changes  of  C0-  IC0-eq  ratio  and  reaction  rate  rB3  are  more  pro¬ 
nounced  at  a  certain  distance  away  from  3PB,  while  at  3PB  (virtually 
within  one-node  distance  to  it)  there  is  almost  no  difference  of  rB3 
between  the  two  cases.  The  B3  reaction  near  x  =  0  is  thus  strongly 
controlled  by  the  boundary  condition  and  relatively  independent 
of  material  property  effects. 

Fig.  12  results  indicate  that  the  increase  of  C0-  eq  can  reduce 
C0-/C0-  eq  while  increasing  rB3.  In  Eq.  (16),  a  reduction  of 
C0-/C0-  eq  and  increase  of  rB3  yield  larger  Cv/Cv,e q  (e.g.  absolute 
Cv)  for  a  given  overpotential.  In  this  way  the  observed  increase 
of  Cy/Cv.eq  at  higher  C0-  eq  together  with  consequent  extension  of 
active  reaction  zone  in  Fig.  1 1  can  be  explained  by  relative  change  of 
C0-/C0-  eq  and  rB3  in  Fig.  12.  An  increase  of  C0-  eq  leads  to  a  higher 
background  chemical  potential,  and  even  though  absolute  concen¬ 
tration  of  surface  oxygen  ion  rises,  the  model  predicts  that  the  local 


Fig.  12.  Distribution  profiles  of  surface  oxygen  concentrations  and  step-B3  reaction  rates  at  -0.2  V  overpotential  in  Case  1  and  Case  2  as  Cv,eq  =  1  x  10-8  mol  cm-3. 
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Fig.  13.  Dependence  of  rates  of  step  S2  and  B3  reactions  at  0.5  p.m  away 
from  3PB  on  overpotential  in  Case  1  (open  symbol)  and  2  (filled  symbol)  as 
Cv.eq  =  1x1 0-8  mol  cm-3 . 


chemical  potential  relating  to  C0-/C0-  eq  drops  (which  is  also  the 
case  for  other  cathodic  overpotentials;  results  not  shown  here).  This 
condition  can  be  expected  to  slow  down  reaction  B3  according  to 
Eq.  (16)  initially,  and  hence  promote  an  increase  of  local  vacancy 
concentration.  Consequently  additional  vacancy  formed  via  2PB 
exchange  reaction  B4  could  be  transported  to  regions  away  from 
2PB  until  new  steady-state  level  of  r 33  is  achieved.  In  this  way  the 
active  reaction  zone  expands,  and  electrode  kinetics  also  improve 
as  C0-  eq  increases. 


4.3.2.  Change  of  surface  reaction  rates  rS2  and  rB3 

As  surface  reaction  rates  rs 2  and  rB 3  determine  local  species 
diffusion  fluxes,  their  profiles  as  functions  of  overpotential  and 
C0-  eq  are  used  to  identify  reason  for  occurrence  of  surface 
diffusion  limit  and  change  of  active  reaction  zone  during  activa¬ 
tion. 

In  Figs.  12  and  13,  rB3  outside  of  3PB  was  more  positive  at  neg¬ 
ative  overpotential,  suggesting  that  the  surface  reaction  is  more 
controlled  by  the  chemical  terms  Cv/Cv,eq  and  C0-/C0-  eq  in  Eq. 
(16)  than  by  overpotential.  The  backward  rate  in  Eq.  (16)  should 
increase  at  more  negative  overpotential.  For  rB3  to  become  more 
positive,  Cv/CVieq  and  C0-/C0-  eq  have  to  increase  accordingly  to 
offset  the  larger  backward  rate.  The  magnified  comparison  in  Fig.  1 4 
shows  that  forward  rates  of  rB3  are  larger  than  backward  rates 
at  x  =  0.5  |jim.  Therefore,  step  B3  reaction  is  confirmed  to  be  more 
“chemically”  driven  outside  of  3PB.  However,  this  result  does  not 
necessarily  mean  that  the  surface  reactions  are  chemical  in  nature. 
Instead  rs 2  remains  2-3  orders  of  magnitude  higher  than  rB3  in 
Fig.  13.  Although  combination  of  reaction  S2  and  B3  would  yield 
zero  charge  transfer  and  hence  a  chemical  oxygen  exchange  reac¬ 
tion  by  formulism,  the  significant  rate  difference  between  S2  and 
B3  reactions  makes  oxygen  intermediates  (O-)  a  stable  phase  dur¬ 
ing  reduction,  and  hence  introducing  the  electrochemical  kinetic 
control  [37]. 

Fig.  13  also  shows  that  at  a  higher  value  of  C0-  eq,  rs 2  increases 
faster  than  rB3,  and  results  in  a  larger  rate  difference  between  rs 2 
and  rB3  for  Case  2  with  C0-  eq  =  1  x  10_10mol  cm-2.  In  the  govern¬ 
ing  equations,  the  increase  of  surface  diffusion  flux  can  be  expected 
from  a  larger  local  rate  difference.  This  explains  the  change  of  3PB 
pathway  kinetics  from  local  reaction  control  for  Case  1  to  surface 
diffusion  control  for  Case  2,  as  shown  by  Figs.  6  and  9.  It  can  be 
also  noticed  in  Fig.  14  that  there  is  an  abrupt  increase  of  rB3  above 
the  overpotentials  for  pathway  transition,  as  is  shown  by  slope 
change  of  the  RB3  plots.  This  indicates  that  abrupt  increase  of  rB3 
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Fig.  14.  Dependence  of  rates  of  backward  and  forward  reaction  rates  of  step  B3  reactions  at  0.5  p,m  away  from  3PB  on  overpotential  in  Case  1  (open  symbol)  and  2  (filled 
symbol)  as  Cv.eq  =  1  x  10-8  mol  cm-3. 
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contributes  to  the  limitation  of  surface  diffusion  and  eventually  the 
3PB-to-2PB  transition. 

4.3.3.  Association  with  cathode  microstructural  design 

Together,  the  above  results  and  discussions  are  consistent  with 
the  assumed  existence  of  an  active  reaction  zone  in  the  MIEC 
cathode  with  micron-level  thickness.  A  change  in  surface  material 
properties  affects  the  3PB-to-2PB  pathway  transition  by  modifying 
the  distribution  of  surface  oxygen  ion  and  bulk  vacancy  and  the  cor¬ 
responding  local  reaction  rates  away  from  the  2PB/3PB  interfaces. 
In  comparison,  Kuznecov  has  observed  that  the  diffusion  length 
of  oxygen  vacancies  from  the  YSZ  electrolyte  into  the  LSM  cath¬ 
ode  was  about  2  f±m  at  950  °C  [49],  while  in  general  the  reported 
effective  oxygen  diffusion  thickness  for  LSM  derived  from  oxygen 
chemical  diffusivity  (Dchem)  and  surface  exchange  coefficient  (K) 
ranges  from  [xm  to  cm-level  [42-44].  This  large  discrepancy  may 
reflect  sophisticated  control  over  data  acquisition  and  processing 
required  for  thermodynamic  property  measurement.  Simulation  in 
this  model  does  predict  a  discernible  extension  of  oxygen  vacancy 
diffusion  depth,  for  improved  surface  catalytic  property  as  a  higher 
C0-  eq  value  implies  more  surface  reaction  sites  available.  This 
finding  may  offer  some  guidance  on  the  cathode  performance 
improvement  methodology  through  microstructure  design.  It  is 
a  conventional  route  to  improve  kinetics  of  a  poor  MIEC  cathode 
through  infiltration  of  ionically  conductive  second-phase  particles 
such  as  doped  ceria.  However,  the  confined  active  reaction  zone 
near  the  electrochemical  interfaces  in  Case  1  implies  that  impreg¬ 
nated  particles  need  to  be  in  good  contact  with  the  3PB  or  extended 
electrolyte  network  out  from  the  2PB  to  effectively  affect  reaction 
rates.  Thus,  control  over  processing  parameters  such  as  precur¬ 
sor  viscosity,  material  loading  and  sintering  temperature  would  be 
important  to  optimize  the  particle  size  and  distribution.  In  contrast, 
the  Case  2  study  indicates  that  impregnation  with  a  highly  cat¬ 
alytic  material  could  improve  surface  reaction  kinetics  and  extend 
the  active  reaction  zone  at  moderate  overpotentials.  Hence  poten¬ 
tially  the  second  methodology  may  benefit  from  less  restriction  on 
secondary  particle  size  and  distribution  to  achieve  similar  perfor¬ 
mance  improvement,  since  contribution  from  surface  diffusion  and 
reactions  outside  3PB  becomes  more  dominant. 

5.  Model  limitation 

The  oxygen  reduction  process  on  the  MIEC  cathode  is  broken 
down  into  six  elementary  steps  (SI  -B4)  for  the  two  pathways  in  this 
model,  and  we  tentatively  propose  step  S2,  B3  and  B4  as  RLSs  based 
on  literature  survey  and  boundary  conditions  required  for  numer¬ 
ical  solution.  Although  these  are  fair  assumptions  at  the  given 
system  conditions  according  to  experimental  results  [10,13,22], 
there  are  other  possible  oxygen  reduction  scenarios  with  different 
elementary  steps  and  RLSs.  A  recent  study  on  LSM  microelectrodes 
by  MIT  research  group  has  assigned  charge  transfer  at  2PB/3PB 
interfaces  as  RLS  below  700  °C  and  electrode  surface  reaction  as 
RLS  above  700  °C  [24].  This  finding  suggests  that  there  is  room 
for  adjustment  of  the  RLS  and  corresponding  boundary  conditions 
at  temperatures  of  interest.  Lately,  quantum  chemistry  and  first 
principle  analysis  show  that  oxygen  adsorption  is  accompanied 
by  formation  of  molecular  oxygen  ions  at  the  electrode  surface 
(O2  and  ol~),  which  can  be  subsequently  converted  into  O- 
[28,36].  This  prediction  would  complicate  RLS  assignments  and 
governing  equations  derivation  by  introducing  additional  elemen¬ 
tary  steps.  Another  important  factor  that  needs  to  be  taken  into 
consideration  for  future  improvement  is,  as  pointed  out  by  Fleig 
and  Mebane  [27,28],  that  the  change  of  surface  potential  across 
the  gas/MIEC  interface  may  in  general  deviate  from  the  nominal 
bulk  overpotential  across  MIEC/electrolyte  interface.  Thus,  it  would 


be  necessary  to  decouple  overpotential  driving  surface  reactions 
from  that  driving  bulk  vacancy  exchange.  Fleig’s  study  showed  that 
surface  overpotential  would  either  linearly  depend  on  the  bulk 
overpotential  or  remain  constant  with  the  range  between  0  and 
2rj  ( rj :  bulk  overpotential)  [27].  In  this  model,  surface  overpoten¬ 
tial  is  assumed  to  be  equal  to  bulk  overpotential  for  computational 
simplification.  To  overcome  this  complexity,  a  detailed  relation¬ 
ship  between  surface/bulk  overpotentials  has  to  be  derived  with 
corresponding  RLS  assumptions  of  surface  reactions. 

6.  Conclusion  and  future  work 

In  this  study  multi-step  charge  transfer  is  incorporated  into  a 
multi-domain  1-D  physical  model  designed  for  the  oxygen  reduc¬ 
tion  scenario  to  examine  the  competitive  behaviors  between  the 
paralleled  3PB  and  2PB  kinetic  pathway  during  MIEC  electrode  acti¬ 
vation.  Analysis  by  V-I  polarization  curve,  Tafel  estimation  and  local 
3PB  current  constitution  has  identified  the  limitation  of  surface 
oxygen  ion  diffusion  as  the  mechanism  for  3PB-to-2PB  kinetic  tran¬ 
sition.  The  transition  voltage  between  -0.2  and  -0.4  V  depends  on 
both  material  surface  and  bulk  parameters.  The  transition  voltage 
becomes  more  negative  at  higher  C0-  eq  and  less  negative  at  higher 
Cv,eq.  The  model  also  demonstrated  surface  reactions  are  driven 
predominantly  by  electrochemical  forces  at  the  3PB,  while  being 
controlled  by  oxygen  vacancy  concentration  variation  at  regions 
away  from  3PB.  Consequently,  improving  of  material  surface  cat¬ 
alytic  activity  can  effectively  promote  electrode  performance  by 
extending  the  active  reaction  zone  outside  of  the  3PB/2PB  interface, 
and  may  be  more  efficient  than  addition  of  a  separate  ionic  con¬ 
ductor  into  the  electrode  via  methods  such  as  infiltration.  Future 
development  of  the  model  should  include:  (1)  adoption  of  a  phys¬ 
ically  more  plausible  oxygen  adsorption  scenario  into  1-D  model 
involving  peroxide  and  superoxide  ions;  (2)  incorporate  into  the 
governing  reaction  mechanisms  the  interrelationship  between  sur¬ 
face  potential  change  and  bulk  overpotential  derived  from  specific 
assumptions  on  the  surface  charging  layer  configuration;  (3)  derive 
cathode  impedance  from  transient  simulation  results  with  defined 
system  properties  of  material,  temperature  and  oxygen  partial 
pressure. 
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Appendix  A.  Parameter  estimation 

Oxygen  vacancy  diffusivity  Dv,  which  is  Dbchem  in  this  model, 
is  reported  by  Mizusaki  to  be  almost  independent  of  perovskite 
compositions  and  vacancy  concentrations  at  700-1 000  °C  [50]. 
Yasuda  found  that  Dv  for  Lao.8Sr0.2Mn03  (LSM20)  with  values  of 
10-6-10-5  cm-2  s-1  at  850  °C  varied  little  with  oxygen  partial  pres¬ 
sure  (p02)[44].  Thus  Dbchem  is  set  as  1  x  10-6  cm-2  s-1  for  LSM  at 
800  °C.  Diffusivity  of  surface  oxygen  ion  Ds  chem  is  kept  the  same 
as  Db  chem  [33,35].  Oxygen  vacancy  concentration,  which  is  sensi¬ 
tive  to  material  composition,  p02  and  temperature,  is  estimated  by 
oxygen  tracer  diffusivity  D*  [42,46]: 

Dq  =/Dv[Vq  ]  (Al) 
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where  the  correlation  factor /is  taken  as  0.69  for  perovskite-type 
crystal.  D*  reported  for  LSM20  at  800°C  ranges  from  4x  10-15 
to  1  x  10-13  for  p02  between  0.13  and  1  atm  [42-45].  Substitu¬ 
tion  of  Dv  and  D*  into  Eq.  (Al)  gives  vacancy  concentration  of 
5.8  x  10-9  to  1.4  x  10-7.  Thus  Cv,MiEC,eq  is  chosen  as  1  x  10-8  to 
1  x  10-7  mol  cm-3,  close  to  Yasuda’s  value  for  LSM20  at  900°C 
[45].  Literature  data  are  rare  for  the  equilibrium  concentration  of 
surface  oxygen  ion,  but  this  value  can  be  estimated  from  the  equi¬ 
librium  oxygen  coverage  and  surface  adsorption  site  density  as  r0 
eq.  r  is  reported  to  be  10-8-10-9  mol  cm-2  for  oxygen  adsorp¬ 
tion  and  0e q  is  usually  assumed  to  be  0.01-0.1  [34,49].  So  C0-  ea 

a<r4 

is  set  as  1  x  10-1°  to  1  x  10-11  mol  cm-2.  Volume-specific  surface 
area  depends  on  average  particle  size  of  the  electrode  according  to 
close-packaging  model  [51]: 


AS 

AV 


6 

davg 


(1  -e) 


(A2) 


where  davg  is  the  average  particle  size,  and  s  is  0.26  for  face- 
centered  cubic  lattice.  Then  AS/  A  V=  8.9  x  1 06  nrr 1  if  daVg  =  0.5  p,m. 
So  AS/AV  is  set  as  1  x  105  cm-1  for  submicron-sized  particles. 

The  exchange  reaction  rates  rs 2,o.  tB3,o  and  rB4,o  directly  relate  to 
the  exchange  current  densities.  AC  impedance  and  DC  polarization 
tests  are  frequently  used  to  determine  i0  for  the  specific  RLS,  how¬ 
ever  the  reported  values  vary  due  to  different  testing  conditions 
and  disagreement  on  kinetic  interpretation,  van  Heuveln  obtained 
i0  of  20-1 10  mA cm-2  at  950  °C  for  a  LSM  electrode  with  surface 
electron  donation  and  surface  oxygen  diffusion  as  RLSs  [10].  While 
the  same  RLSs  were  also  assigned  by  S.  Wang  on  LSM,  i0  was  set  as 
1 0  mA  cm-2  at  900  °C  [  1 3  ].  Chen’s  fitting  of  i0  yielded  1-10  mA  cm-2 
for  LSM  under  much  lower  oxygen  partial  pressure  at  800  °C  [48]. 
At  even  lower  temperature  of  700  °C  Horita  still  obtained  an  i0  of 
2  mAcm-2  for  a  LSM  mesh  electrode  on  YSZ  [23].  The  value  of  rS2,o 
selected  in  this  model  gives  z'o  around  10 mAcm-2.  Estimation  of 
rB3,o  and  rB4,o  is  difficult  due  to  lack  of  data  for  a  LSM  electrode. 
The  corresponding  rates  in  Svensson  and  Coffey’s  models  convert 
toi0  well  below  1  mAcm-2  [34,35].  Such  a  low  value  means  that  the 
enhancement  of  electrode  kinetics  from  2PB  pathway  may  be  rather 
limited.  However,  if  LSM  does  become  a  good  MIEC  like  LSF  and  LSC 
as  expected  at  higher  overpotential,  it  should  exhibit  similar  kinetic 
properties  for  oxygen  reduction  as  these  materials.  The  exchange 
current  densities  assigned  for  LSC  and  LSCF  are  generally  between 
330  and  1100  mA  cm-2  at  800  °C  [52,53  ].  Hence,  rB30  is  set  equal  to 
rS2,o  considering  that  surface  exchange  coefficiency  of  LSM  is  about 
1-2  orders  lower  than  other  good  MIECs,  and  essentially  3PB  and 
2PB  pathways  would  be  competitive  on  surface  reactions.  The  value 
of  rB40  is  chosen  so  that  iB4)0  equal  to  193  mAcm-2,  indicating  step 
B4  can  be  RLS  at  more  negative  overpotential. 


Appendix  B.  Potential  derivation 

In  Svensson  and  Coffey’s  models,  the  driving  force  for  2PB  and 
3PB  reactions  are  defined  as  the  electrochemical  potential  differ¬ 
ence  between  the  involved  reaction  species  across  the  interface, 
which  for  2PB  reaction  can  be  expressed  as  [34,35]: 


9FV,S  _  ..MIEC  ..YSZ 

zr//2PB  “  ^v,2PB  ^v,2PB 


(Bl) 


where  Z78PB  is  the  defined  driving  force  for  2PB  reaction,  and 
and  /jtf s2pB  are  the  electrochemical  potentials  of  oxygen  vacancy 
in  MIEC  cathode  and  electrolyte,  respectively.  Then  the  measured 
total  potential  drop  between  WE  and  RE  in  the  3-electrode  system 


is  defined  in  the  models  from  2PB  potential  deviation  from  OCV  as 
[34,35]: 


rMIEC 

S,2PB 


p  _  pOC  RT  Ly,2PB  1  /  MIEC  ..YSZ  \ 

L2pb  -L2PB  2F  in  MIEC  +  2P^v,2PB  ^v,2PB' 
Lv,2PB 


(B2) 


where  C°c^BIEC  and  are  the  2PB  oxygen  vacancy  concentration 

in  MIEC  cathode  at  OCV  and  polarization,  respectively.  For  a  1-D 
model  assuming  uniform  lateral  potential  distribution,  the  mea¬ 
sured  3PB  and  2PB  potentials  must  be  the  same,  which  is  also  the 
case  at  OCV,  hence  [35]: 


c  _  r  r  nOC  n  rOC 

L2PB  —  L3PE},  L2pB  ^2PB  —  ^3PB  ^3PB 


(B3) 


In  our  model,  we  define  this  measured  potential  difference  from 
OCV  as  the  “net”  driving  forces  for  both  3PB  and  2PB  reactions, 
therefore  we  have: 


??2PB  =  ??3PB  (B4) 

The  physical  implication  of  this  overpotential  definition  is 
the  deviation  of  the  electrostatic  potential  difference  across  the 
interfaces  from  OCV.  With  the  assumption  of  dilute  solution,  the 
electrochemical  potentials  for  charged  species  in  Eq.  (B2)  can  be 
expressed  as  [34]: 

..MIEC  //0  ,  nT  in  rMIEC  ,  nnvfcMIEC  ..MIEC  ,  nnvfcMIEC 

^i/,2PB  =  ln  S,,2PB  +  Ztq>2PB  =  ^y,chem,2PB  +  Zbq>2PB 

(B5) 


/4?PB  =  >4  +  RT  ln  CJS2PB  +  2F0lm  =  <Schem,2PB  +  2F0lm  (BB) 

where  fi„  is  the  stand  reference  state,  0  the  electrostatic  poten¬ 
tial  and  fiv  chem  the  chemical  potential  of  oxygen  vacancy  at  2PB. 
Substitute  Eqs.  (B5)  and  (B6)  into  Eq.  (B2),  it  can  written  that 


c  cOC  RT  , 

e2pb  =  £2pB  -  2P  ln 


HVIIEC  1 
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(B7) 


At  OCV  we  know  that  /x°c ^BEC  =  /x°2pbZ’  which  leads  to 

oc.MIEC  .  9  p(  /jsoc.MIEC  _  ^oc.YSZn  _  oc.YSZ 
^V,chem,2PB  +  zr^2PB  ^2PB  '  ~  ^y,chem,2PB 


(B8) 


And  if  vacancy  concentration  in  electrolyte  is  assumed  to  keep 
constant  during  reaction,  then  we  get  /x^lem  2PB  =  4S.2PB-  So 
Eq.  (B8)  can  be  substituted  into  Eq.  (B7)  to  get 


E2PB  =  E2pb  -  2F  In 


rMIEC  . 
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(B9) 


Therefore,  by  eliminating  the  vacancy  concentration  terms  from 
Eq.  (B9)  we  obtain 

>)2PB  =  E2PB  -  E2PB  =  «PBC  ~  ^IpD  ~  (B1«) 

With  introduction  of  surface  oxygen  intermediate  ion  O-,  we 
can  also  derive  the  modified  driving  forces  for  the  3PB  reaction  as: 

?73PB  =  +  ^3PMB1EC  -  2^S)  -  (^PBMEC  +  ^3°PB  MIEC  -  2^3^)  (B1  1  ) 

where  is  the  surface  electrostatic  potential  of  MIEC  at  3PB. 

If  zero  electrical  field  is  assumed  between  3PB  and  2PB  interfaces, 
then  we  have 

^BK  =  «C.  <*>3TO  =  0lm  (B12) 
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The  derivation  of  r\2 pb  =  *13pb  would  imply  03,^IEC  = 
which  indicates  that  at  triple-phase-boundary  where  the  MIEC 
cathode  surface  emerges  with  electrolyte,  its  electrostatic  poten¬ 
tial  would  assume  the  same  value  as  that  of  electrolyte  surface. 
With  similar  definition  from  Eqs.  (BIO)  and  (Bll),  we  can  obtain 
the  surface  overpotential  r\s  as 

r]s=Es-  £°C  =  (0MIEC  -  0S-MIEC)  _  (0OC.MIEC  _  ^oc.S.MIEC  j  (B13) 

where  surface  electrode  potential  Es  is  difference  between  the  bulk 
electrostatic  potential  of  MIEC  <£>MIEC  and  the  surface  electrostatic 
potential  of  MIEC  <£S-MIEC.  The  change  of  electrostatic  potential  dif¬ 
ference  in  Eqs.  (BIO),  (Bll)  and  (B13)  as  net  driving  forces  for 
electrochemical  reactions  can  then  be  used  to  obtain  B-V  kinetic 
Eqs.  ( 1 4)-(  1 7)  through  transition-state-theory  analysis  in  Mebane’s 
model  [28].  As  discussed  in  Sections  2.3  and  5  the  model  assumed  a 
simplified  relationship  between  surface  potential  step  and  interfa¬ 
cial  overpotentials,  out  from  zero  electrical  field  assumption  inside 
the  MIEC  bulk  and  on  its  surface.  For  a  poor  ionic  MIEC  like  LSM,  its 
bulk  background  charge  would  be  5-6  orders  of  that  carried  by  oxy¬ 
gen  vacancy  with  initial  comparison  of  this  study  to  Mebane  and 
Svensson’s  models  [34,37].  For  the  simulated  overpotential  range 
in  this  model,  the  increase  of  bulk  charge  due  to  vacancy  exchange 
into  MIEC  can  be  estimated  to  be  around  1%  at  most.  Therefore  as 
estimated  by  Svensson  et  al.  [34],  the  species  migration  flux  in  this 
case  can  be  negligible  to  the  diffusion  flux.  However,  the  model  can 
break  down  for  simulation  of  good  ionic  MIEC  with  much  higher 
intrinsic  oxygen  vacancy  concentration.  And  then  a  modified  model 
with  independent  consideration  of  surface  potential  and  electrical 
field  effects  would  be  required. 
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